import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

# R(t) parameters
a_r = 3.7741
b_r = 0.5769
c_r = -0.8062
d_r = 0.0811

# A(t) parameters
a_a = 0.6509
b_a = -0.2017
c_a = -1.1990
d_a = -0.4092

# M(t) parameters
a_m = 0.120563
b_m = 0.546877
c_m = -242.873714

# V(t) parameters
a_v = -346.3630
b_v = -0.1895
c_v = -0.9567
d_v = 0.1720

# W(N_t, C, k) parameters
alpha = 3
lamda = 1.5
mu = 1.5
k = 0.05

visitor_mean =  991194.4444444445
visitor_std =  286175.79770371574
satisfaction_mean =  0.3744444444444444
satisfaction_std =  0.020631990023884737
expanditure_mean =  261884.22222222222
expanditure_std =  62927.576153115675
tax_mean =  1447762.9444444445
tax_std =  549045.7156859094
investment_mean =  1763810.0555555555
investment_std =  714385.1499613319
glacier_mean =  4.331388888888889
glacier_std =  1.1985590950474607

# 安全的指数函数
def safe_exp(x):
    return np.exp(np.clip(x, -10, 10))  # 限制指数函数的输入范围

def R(N_t, C):
    result = N_t * a_r * safe_exp(c_r + b_r * C) - d_r * C
    return result

def F(N_t):
    result = N_t * 16 / 365
    return result

def V(N_t, C, t):
    # 对 N_t 进行适当的缩放以避免溢出
    N_t_scaled = N_t
    result = a_v + b_v * C * safe_exp(c_v * N_t_scaled) + d_v * t
    return result

def A(N_t, C, t):
    V_value = V(N_t, C, t)
    result = a_a + b_a * safe_exp(C) + c_a * N_t + d_a * V_value
    return result

def M(N_t, t):
    result = c_m + b_m * N_t + a_m * t
    return result

def W(N_t, C, k, t, alpha, lamda, mu):
    R_value = R(N_t, C)
    F_value = F(N_t)
    A_value = A(N_t, C, t)
    M_value = M(N_t, t)
    result = alpha * R_value - lamda * F_value + mu * A_value + k * alpha * M_value
    return result

# 读取数据
file_path = './assets/normalized_data(2).xlsx'  # Excel 文件路径
data = pd.read_excel(file_path)

# 提取年份数据
years = data['T'].values
N_values = data['N1'].values
C_values = data['C1'].values

# 计算原本的 W 曲线
original_results = []
for t, N, C in zip(years, N_values, C_values):
    W_original = W(N, C, k, t, alpha, lamda, mu)
    original_results.append(W_original)

# 定义参数变化范围
change_rates = np.linspace(0.8, 1.2, 5)  # 例如，变化范围为 ±20%

# 存储灵敏度分析结果
N_results = []
C_results = []
k_results = []

# 计算不同参数值下的 W 曲线
for change in change_rates:
    N_results.append([W(N * change, C, k, t, alpha, lamda, mu) for t, N, C in zip(years, N_values, C_values)])
    C_results.append([W(N, C * change, k, t, alpha, lamda, mu) for t, N, C in zip(years, N_values, C_values)])
    k_results.append([W(N, C, k * change, t, alpha, lamda, mu) for t, N, C in zip(years, N_values, C_values)])

# 绘制灵敏度分析图
sns.set_theme(style="darkgrid")
palette = sns.color_palette("PuBu", 5)

# 绘制 N_t 变化的 W 曲线
plt.figure(figsize=(10, 6))
sns.lineplot(x=years, y=original_results, label='Original W', color='black')
for i, change in enumerate(change_rates):
    sns.lineplot(x=years, y=N_results[i], label=f'N_t {int((change-1)*100)}%', color=palette[i])
plt.ylabel('Total Profit', fontsize=14, fontname='Times New Roman')
plt.title('Sensitivity Analysis of N_t', fontsize=16, fontname='Times New Roman')
plt.legend()
plt.xticks(ticks=np.arange(years.min(), years.max() + 1, 1), fontname='Times New Roman')  # 设置横轴刻度为整数
plt.yticks(fontname='Times New Roman')
plt.show()

# 绘制 C 变化的 W 曲线
plt.figure(figsize=(10, 6))
sns.lineplot(x=years, y=original_results, label='Original W', color='black')
for i, change in enumerate(change_rates):
    sns.lineplot(x=years, y=C_results[i], label=f'C {int((change-1)*100)}%', color=palette[i])
plt.ylabel('Total Profit', fontsize=14, fontname='Times New Roman')
plt.title('Sensitivity Analysis of C', fontsize=16, fontname='Times New Roman')
plt.legend()
plt.xticks(ticks=np.arange(years.min(), years.max() + 1, 1), fontname='Times New Roman')  # 设置横轴刻度为整数
plt.yticks(fontname='Times New Roman')
plt.show()

# 绘制 k 变化的 W 曲线
plt.figure(figsize=(10, 6))
sns.lineplot(x=years, y=original_results, label='Original W', color='black')
for i, change in enumerate(change_rates):
    sns.lineplot(x=years, y=k_results[i], label=f'k {int((change-1)*100)}%', color=palette[i])
plt.ylabel('Total Profit', fontsize=14, fontname='Times New Roman')
plt.title('Sensitivity Analysis of k', fontsize=16, fontname='Times New Roman')
plt.legend()
plt.xticks(ticks=np.arange(years.min(), years.max() + 1, 1), fontname='Times New Roman')  # 设置横轴刻度为整数
plt.yticks(fontname='Times New Roman')
plt.show()